; void ondas_asm (
; 	unsigned char *src, RDI
; 	unsigned char *dst, RSI
; 	int m,				EDX
; 	int n,				ECX
; 	int row_size		r8d
;	int y0				r9d
;	int x0				[RSP+8] o despues del stackFrame es [rbp+16]
; );

global ondas_asm

section .rodata

PI: DD 3.1415,3.1415,3.1415,3.1415
RADIUS: DD 35.0,35.0,35.0,35.0
WAVELENGTH: DD 64.0,64.0,64.0,64.0
TRAINWIDTH: DD 3.4,3.4,3.4,3.4

UNO: DD 1,1,1,1
UNOF: DD 1.0, 1.0, 1.0, 1.0
DOS: DD 2.0,2.0,2.0,2.0
SEIS: DD 6.0,6.0,6.0,6.0
SEISCUATRO: DD 64.0,64.0,64.0,64.0
CIENTOVEINTE: DD 120.0, 120.0, 120.0, 120.0
CINCUENTACUARENTA: DD 5040.0,5040.0,5040.0,5040.0


section .text

ondas_asm:
	push rbp
	mov rbp, rsp
	push rbx
	push r12
	push r13
	push r14
	push r15
		
	pxor xmm14, xmm14		;Limpiamos xmm14
	movd xmm14, r9d			;xmm14 = 0|0|0|x0
	movdqu xmm13, xmm14
	pslldq xmm13, 4			;xmm13 = 0|0|y0|0
	paddd xmm14, xmm13		;xmm14 = 0|0|y0|y0
	movdqu xmm13, xmm14
	pslldq xmm13, 8			;xmm13 = y0|y0|0|0
	paddd xmm14, xmm13		;xmm14 = y0|y0|y0|y0
	
	;Hacemos lo mismo con x0
	pxor xmm15, xmm15		;Limpiamos xmm15
	movd xmm15, [rbp+16]	;xmm15 = 0|0|0|x0
	movdqu xmm13, xmm15
	pslldq xmm13, 4			;xmm13 = 0|0|x0|0
	paddd xmm15, xmm13		;xmm15 = 0|0|x0|x0
	movdqu xmm13, xmm15
	pslldq xmm13, 8			;xmm13 = y0|y0|0|0
	paddd xmm15, xmm13		;xmm15 = x0|x0|x0|x0	
	
	;En xmm14 tenemos x0 repetidos 4 veces y en xmm15 tenemos y0 repetido 4 veces
	
	xor r12, r12 	;Contador filas
	
.cicloFilas:	
	xor r13, r13	;Contador columnas
	
.cicloColumnas:
	mov rbx, r12 	; j
	imul ebx, r8d	; j*row_size
	add rbx, r13	; avanzamos de a 8
	
	;calculamos los dx : calculamos los x para los 8 pixeles
	
	movdqu xmm13, [UNO] ;guardamos el uno en xmm13
	
	pxor xmm1, xmm1	;limpiamos xmm1
	paddd xmm1, xmm13 ;sumamos 1 a todo
	pslldq xmm1, 4	;xmm1 = 1|1|1|0
	paddd xmm1, xmm13 ;sumamos 1 a todo
	pslldq xmm1, 4	;xmm1 = 2|2|1|0
	paddd xmm1, xmm13 ;sumamos 1 a todo
	pslldq xmm1, 4	;xmm1 = 3|2|1|0
	
	pxor xmm2, xmm2		;Limpiamos xmm2
	mov r10d, 4			; 4 como entero
	movd xmm2, r10d		;xmm2 = 0|0|0|4
	movdqu xmm3, xmm2
	pslldq xmm3, 4		;xmm3 = 0|0|4|0
	paddd xmm2, xmm3	;xmm2 = 0|0|4|4
	movdqu xmm3, xmm2
	pslldq xmm3, 8		;xmm3 = 4|4|0|0
	paddd xmm2, xmm3	;xmm2 = 4|4|4|4
	
	paddd xmm2, xmm1	;xmm2 = 7|6|5|4
	
	pxor xmm3, xmm3 ;limpiamos xmm3
	movd xmm3, r13d ;xmm3 = 0|0|0|x
	movdqu xmm4, xmm3
	pslldq xmm4, 4	;xmm4 = 0|0|x|0
	paddd xmm3, xmm4 ;xmm3 = 0|0|x|x
	movdqu xmm4, xmm3
	pslldq xmm4, 8	;xmm4 = x|x|0|0
	paddd xmm3, xmm4 ;xmm3 = x|x|x|x
	
	paddd xmm1, xmm3	;xmm1 = x+3|x+2|x+1|x
	paddd xmm2, xmm3	;xmm2 = x+7|x+6|x+5|x+4
	
	;xmm1 y xmm2 tienen los x
	
	psubd xmm1, xmm15	;xmm1 = x+3-x0|x+2-x0|x+1-x0|x-x0
	psubd xmm2, xmm15	;xmm2 = x+7-x0|x+6-x0|x+5-x0|x+4-x0
	
;xmm1 y xmm2 tienen los dx
	
	;calculamos los dy : como todos estan en la misma fila todos tienen el mismo y
	
	pxor xmm3, xmm3 ;limpiamos xmm5
	movd xmm3, r12d ;xmm3 = 0|0|0|y
	movdqu xmm4, xmm3
	pslldq xmm4, 4	;xmm4 = 0|0|y|0
	paddd xmm3, xmm4 ;xmm3 = 0|0|y|y
	movdqu xmm4, xmm3
	pslldq xmm4, 8	;xmm4 = y|y|0|0
	paddd xmm3, xmm4 ;xmm3 = y|y|y|y
	
	;xmm3 tiene el y
	
	psubd xmm3, xmm14	;xmm3 = y-y0|y-y0|y-y0|y-y0
	
;xmm3 tiene los dy
	
	cvtdq2ps xmm1, xmm1
	cvtdq2ps xmm2, xmm2 ; convertimos a float los dx
	
	cvtdq2ps xmm3, xmm3 ; convertimos a float los dy
	
	mulps xmm1, xmm1
	mulps xmm2, xmm2 ; elevamos al cuadrado dx
	
	mulps xmm3, xmm3 ; elevamos al cuadrado dy
	
	addps xmm1, xmm3
	addps xmm2, xmm3 ; sumamos los cuadrados dx*dx + dy*dy
	
	sqrtps xmm1, xmm1
	sqrtps xmm2, xmm2 ;hacemos la raiz cuadrada sqrt(dx*dx + dy*dy)
	
;xmm1 y xmm2 guardan dxy como floats
	
	movdqu xmm3, [RADIUS] ;guardamos el radio en xmm3
	movdqu xmm4, [WAVELENGTH] ;guardamos el wavelength en xmm4
	
	subps xmm1, xmm3
	subps xmm2, xmm3 ;restamos el radio a dxy
	
	divps xmm1, xmm4
	divps xmm2, xmm4 ;dividimos (dxy-radius) por el wavelength
	
	movdqu xmm3, xmm1
	movdqu xmm4, xmm2 ; guardamos una copia de r
	
;xmm1 y xmm2 guardan r, y xmm3 y xmm4 guardan una copia
	
	roundps xmm5, xmm1, 01b
	roundps xmm6, xmm2, 01b ; calculamos floor(r), 01b redondea hacia abajo
	
;xmm5 y xmm6 guardan floor(r)
	
	movdqu xmm7, [TRAINWIDTH] ;guardamos el trainwidth en xmm7
	
	divps xmm1, xmm7
	divps xmm2, xmm7 ;dividimos r por el trainwidth
		
	mulps xmm1, xmm1
	mulps xmm2, xmm2 ;elevamos al cuadrado
	
	movdqu xmm13, [UNO] ;guardamos el uno en xmm13
	cvtdq2ps xmm13, xmm13	;convertimos a float
	
	addps xmm1, xmm13
	addps xmm2, xmm13 ;sumamos 1
	
	movdqu xmm8, xmm13
	movdqu xmm9, xmm13 ;nos copiamos el uno
	
	divps xmm8, xmm1
	divps xmm9, xmm2 ;calculamos la inversa a
	
	movdqu xmm1, xmm8
	movdqu xmm2, xmm9 ; pasamos a xmm1 y xmm2 la inversa a
	
;xmm1 y xmm2 guardan la inversa a
	
	subps xmm3, xmm5
	subps xmm4, xmm6 ; restamos floor(r) a r
	
	movdqu xmm5, [DOS] ; guardamos el dos en xmm5
	
	movdqu xmm6, [PI] ;guardamos PI en xmm5
	
	mulps xmm5, xmm6 ; multiplicamos 2*PI
	
	mulps xmm3, xmm5
	mulps xmm4, xmm5 ; multiplicamos (r-floor(r))*2*PI
	
	subps xmm3, xmm6
	subps xmm4, xmm6 ; restamos PI
	
;xmm3 y xmm4 guardan t

	movdqu xmm5, xmm3
	movdqu xmm6, xmm4 ; guardamos una copia de t
	
	mulps xmm5, xmm3
	mulps xmm6, xmm4 ; calculamos t*t
	
	mulps xmm5, xmm3
	mulps xmm6, xmm4 ; calculamos t*t*t
	
	movdqu xmm11, xmm5
	movdqu xmm12, xmm6 ;Tengo una copia de t*t*t
	
	movdqu xmm13, [SEIS]
	
	divps xmm5, xmm13
	divps xmm6, xmm13	;t*t*t* / 6
	
	mulps xmm11, xmm3
	mulps xmm12, xmm4	;t*t*t*t
	
	mulps xmm11, xmm3
	mulps xmm12, xmm4	;t*t*t*t*t 
	
	movdqu xmm7, xmm11
	movdqu xmm8, xmm12  ;t*t*t*t*t
	
	movdqu xmm13, [CIENTOVEINTE]
	
	divps xmm7, xmm13
	divps xmm8, xmm13 ;t*t*t*t*t / 120
	
	mulps xmm11, xmm3
	mulps xmm12, xmm4	;t*t*t*t*t*t
	
	mulps xmm11, xmm3
	mulps xmm12, xmm4	;t*t*t*t*t*t*t
	
	movdqu xmm9, xmm11
	movdqu xmm10, xmm12 ;t*t*t*t*t*t*t
	
	movdqu xmm13, [CINCUENTACUARENTA]
	
	divps xmm9, xmm13
	divps xmm10, xmm13 ;t*t*t*t*t*t*t / 5040
	
	subps xmm3, xmm5
	subps xmm4, xmm6 ; t - t*t*t/6
	
	addps xmm3, xmm7
	addps xmm4, xmm8 ; t - t*t*t/6 + t*t*t*t*t/120
	
	subps xmm3, xmm9
	subps xmm4, xmm10; t - t*t*t/6 + t*t*t*t*t/120 - t*t*t*t*t*t*t/5040
	
	mulps xmm1, xmm3
	mulps xmm2, xmm4 ; prof = a . ( t - t*t*t/6 + t*t*t*t*t/120 - t*t*t*t*t*t*t/5040)
	
;xmm1 y xmm2 guardan prof

	movdqu xmm3, [SEISCUATRO] ; leemos el 64.
	
	mulps xmm1, xmm3
	mulps xmm2, xmm3 ; calculamos prof*64

	pxor xmm3, xmm3
	movq xmm3, [rdi+rbx] ;Agarramos 8 pixeles de la imagen fuente
	
	movdqu xmm4, xmm3 ; copia de xmm3
	psrldq xmm4, 4 ; ponemos los segundos 4 bytes en la parte mas baja
	pmovzxbd xmm4, xmm4 ; convertimos a DW los segundos 4 bytes extendiendo con cero
	
	pmovzxbd xmm3, xmm3 ; convertimos a DW los primeros 4 bytes extendiendo con cero
	
	cvtdq2ps xmm3, xmm3
	cvtdq2ps xmm4, xmm4 ; convertimos a float
	
	addps xmm1, xmm3
	addps xmm2, xmm4 ; calculamos prof*64 + Isrc(x,y)
	
	cvtps2dq xmm1, xmm1
	cvtps2dq xmm2, xmm2 ; convertimos los float a DW
	
	packssdw xmm1, xmm2 ; convertimos a word saturando los negativos a 0
	packuswb xmm1, xmm1 ; convertimos a bytes saturando los numeros mayores a 255 a 255
	
	movq [rsi+rbx], xmm1 ; guardamos los 8 pixeles en la imagen destino
	
.avanzarColumnas:
	add r13, 8 ; avanzamos a las proximas 8 columnas
	cmp r13d, ecx
	je .avanzarFila
.cuantoAtras:
	lea r10, [r13+8]
	cmp r10d, ecx	;si es menor o igual nos quedan por lo menos 8 columnas por copiar
	jle .cicloColumnas
	sub r13, 1
	jmp .cuantoAtras	;Con el ciclo calculamos cuanto tenemos que ir para atras para agarrar justo 8
	
.avanzarFila:
	add r12, 1	; vamos a la siguiente fila
	cmp r12d, edx
	jl .cicloFilas ; nos quedan filas por recortar
	
.fin:
	pop r15
	pop r14
	pop r13
	pop r12
	pop rbx
	pop rbp
	ret
